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The recently proposed discrete unified gas kinetic scheme (DUGKS) is a finite volume method 
for deterministic solution of the Boltzmann model equation with asymptotic preserving property. 

In DUGKS, the numerical flux of the distribution function is determined from a local numerical 
solution of the Boltzmann model equation using an unsplitting approach. The time step and mesh 
resolution are not restricted by the molecular collision time and mean free path. To demonstrate the 
capacity of DUGKS in practical problems, this paper extends the DUGKS to arbitrary unstructured 
meshes. Several tests of both internal and external flows are performed, which include the cavity 
flow ranging from continuum to free molecular regimes, a multiscale flow between two connected 
cavities with a pressure ratio of 10"^, and a high speed flow over a cylinder in slip and transitional 
regimes. The numerical results demonstrate the effectiveness of the DUGKS in simulating multiscale 
flow problems. 
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I. INTRODUCTION 

Gas flows can be classified into different regimes based on the Knudsen number (Kn), which is defined as the ratio 
of the mean free path of the gas to the physical characteristic length. For flows with Kn > 0.001, non-equilibrium 
effects become important and the classical Navier-Stokes-Fourier (NSF) equations fail to describe such effects [T]. The 
Boltzmann equation can serve as a fundamental equation that is valid for the whole range of Knudsen number. 

There are mainly two types of numerical approaches for solving the Boltzmann equation. The first one is the 
widely used direct simulation Monte Garlo (DSMG) method [T], which is the prevailing technique for simulation of 
high-speed rarefied gas flows. However, in DSMG the particle transport and collision processes are decoupled, and the 
cell size and time step are required to be smaller than the mean free path and the particle collision time, respectively. 
For flows in near continuum or continuum flow regime, this requirement will lead to enormous computational costs. 
Another undesired feature of the DSMG is the statistical noise that must be reduced with time consuming sampling 
and averaging, which will lead to even larger computational costs for low speed flows and transient problems [T] . It is 
noted that some efforts have been devoted to reduce statistical noise of the DSMG methods ElE]. The second approach 
is to solve the Boltzmann equation directly using deterministic numerical schemes. The most popular one of this type 
is the Discrete Velocity Method (DVM) or Discrete Ordinate Method (DOM) [IHS], in which the velocity space is 
discretized into a finite set of discrete velocities and the same operator splitting technique as DSMG is employed to 
solve the discrete-velocity kinetic equation HEIII]. Therefore, these methods face the same constraints on time-step 
and cell-size as the DSMG for the continuum and near-continuum flows. Recently, some asymptotic preserving (AP) 
schemes have been proposed in order to overcome these disadvantages {e.g., IHlin]). These schemes have been shown 
to be able to recover the Euler solutions in the continuum limit, but it is still not clear whether the Navier-Stokes 
solutions can be obtained. 

Recently, a unified gas kinetic scheme (UGKS) in finite-volume formulation was constructed for all Knudsen number 
flows [lUHS]- Unlike the traditional DOM or DVM, the particle transport and collision are considered simultaneously 
in UGKS in the update of the discrete distribution function. Consequently, the restriction on the cell size and time 
step is avoided. Therefore, UGKS can be used to simulate entire Knudsen number flows efficiently [14]. 

A novel discrete unified gas kinetic scheme (DUGKS) for multi-regime flows was proposed recently [TSIIIS], which 
shares the same modeling mechanism as the original UGKS m- The main difference lies in the reconstruction of 
the discrete distribution function at cell-interface. In UGKS, the time-dependent interface distribution function is 
determined from the local integral solution of the kinetic equation, while in DUGKS, the distribution function at 
the half time step is determined from a characteristic solution of kinetic equation. This reconstruction includes the 
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coupled effects of particle transport and collision, and makes the updating rule much simplified in comparison with 
the UGKS. 

In previous works [nmi], the DUGKS has been applied to both low speed and high speed non-equilibrium flows 
based on structured meshes. To further demonstrate the potential applications of the DUGKS for practical problems 
involving complex geometries or large flow gradients, in this work we will extend the DUGKS to arbitrary unstructured 
meshes. 

The rest of the paper is organized as following. In Sec. 2 the general procedure of the DUGKS on unstructured 
meshes is presented. In Sec. 3, several numerical examples, including the micro cavity flow, an expansion flow between 
two connected cavities, and the rarefied gas flow past a circular cylinder, are provided to demonstrate the applicability 
of the current method in simulating flows at or covering different regimes. A brief summary is given in the last section. 


II. DISCRETE UNIFIED GAS KINETIC SCHEME 
A. Shakhov model 


The DUGKS is based on the Boltzmann model equation in which the collision operator is approximated by the 
Shakhov model m for monatomic gases. In D dimensional space, the model equation is 

f+|-V/ = -i[/-/^], (1) 

where / = /(t, rj, x) is the velocity distribution function of particles with velocity ^ = (^i,..., in D dimensional 
velocity space at position x = (xi,... ,xd) and time t. Here t) = (^d+i, ... ,^ 3 ) is a vector in a space with dimen¬ 
sionality L = 3 — D, consisting of the rest components of the particle velocity (^ 1 ,^ 2 , ^ 3 ) in 3-dimensional space; 
is the Shakhov equilibrium distribution function given by the Maxwellian distribution function /®'^, plus a heat flux 
correction term as 


f = r 


1 + (1 - Pr) 


eg 

bpRT 


RT 


-5) 


= r+fPr, 


( 2 ) 


where Pr is the Prandtl number and c — ^ — U is the peculiar velocity with U being the macroscopic fluid velocity; 
g is the heat flux vector, R is the specific gas constant, and T is the temperature. The collision time t in Eq. Q is 
related to the dynamic viscosity p and pressure p by r = p/p. The Maxwellian distribution function is given by 

where p is the macroscopic gas density. The conservative macroscopic flow variables W = (p, pU, pE)^ are calculated 
as velocity moments of the distribution function. 


W = J 


(4) 


where -ip = (l,^, -I- and pE = 5 ( 0 ^ + g^) + CyT is the total energy with Cy being the heat capacity at 

constant volume. The heat flux g is defined by 

g = ^ y c(c2 -k g^)/d|d77. (5) 

The distribution function / depends only on £ in D dimensional velocity space and is irrelevant to r). To remove 
the dependence on 77 , two reduced distribution functions can be introduced [1] 


( 6 a) 

( 6 b) 


The macroscopic variables can be computed from these reduced distribution function as 

P = J gdi, pU = J ^gdt P^ = \J 


g{x,i,t) = J f{$,T],x,t)dr], 
h{x,i,t)=j f{i,-n,x,t)dri. 


( 7 ) 
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and the heat flux can be computed as 


9 = ^ y c{c^g + h)d^. (8) 

The evolution equations for the reduced distribution functions can be deduced from Eq. Q as 

^+^. [g-/] , (9a) 

^+^.Vh=ng = -^[h-h^], (9b) 

where the reduced equilibrium distribution functions and can be deduced from the original equilibrium distri¬ 
bution as 


with 


g^{x, J /'^(^ > V, X, t)dr] = -h gpr, (10a) 

h^{x, ^,t) = J T], X, t)dr] = -h hpr, (10b) 
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B. Discrete unified gas kinetic scheme on unstructured meshes 


1. Updating of the cell-averaged distribution function 


The updating rules for g and h in Eq. (101 have the same structure 




( 12 ) 


where (j) = g or h. The generic symbol (j) will be used to denote g and h in the following. The DUGKS is an explicit 
finite volume scheme for solving the kinetic equation (12|. The computation domain is firstly divided into some 


control volumes (cells). By integrating Eq. (12) in each cell from time to t„+i, we have 




(13) 


Here (j)j and Ulj is the cell averaged (j) and fl in cell j; \Vj\ is the cell volume and At = tn+i — tn is the time step. Note 
that the trapezoidal and middle-point rules are used for the collision and convection terms in Eq. (13), respectively. 
The term _ 7 r"+i /2 jg ^ across the interface of cell j and is evaluated as 




(14) 


where Sj is the outward normal vector of the fcth face of cell j with face area and Xj is the center of the face. 
Equation( |13[ ) can be rewritten in an explicit form by introducing the transformed distribution functions nulls], (j) 
and (jy^ 


t>r 






1 - 1 - 1/2 


(15) 
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where 


= (t>+ = 


2t + At 


2t 


2t- At ~ 
2t + At^ 



2At 

2t +At 


(16a) 

(16b) 


Due to the conservative property of the collision term, the conservative variables can also be calculated from the 
transformed distribution functions (j) as m 

P = J ffdl, pU^ J |gd|, pE=\j {e~9 + ^)d|, (17) 

and 

^ ^ 2T+^AfPr ^’ ^ ^ ^ ^ 


Therefore, in actual implementation, the evolution of transformed distribution functions (j) is tracked according to 
Eq. ( |15| ), instead of the original distribution functions (j) in order to avoid implicit computations. This is one of the 
major differences between the DUGKS and the UGKS methods. 


2. Flux evaluation on unstructured mesh 


To update according to Eq. (IT^, the flux T. 


1 + 1/2 . 


is required. From the definition of T, 


■n+l/2 


given by Eq. (141, 


the original distribution functions at middle time step at cell interfaces, ie., ^) is needed. This is done by 

solving the kinetic equation ([T^ locally around the cell interface. To this end, Eq. (12) is integrated from time to 
^n+i /2 along a characteristic line which ends at the face center aij. 


- rixf - Is,I) = - I) + n^ixf - Is, |) 


(19) 


where s = t„+i /2 — tn is the half time step. Here the trapezoidal rule is used again for the collision term. Similar to 
the treatment of Eq. (13), another two transformed distribution functions are introduced as 


s 2r + s 2 a 

> = (h -H =- d> -d)*, 

^2 2t ^ 2t ' 

, s^ 2r — s - 2s 

2 2r + s^ 2r + s ' 


then Eq. (19) can be expressed explicitly as 

,^"+i/2(a;/,|) = <^+'"(a;/-|s,|). 


(20a) 

(20b) 

( 21 ) 


Piecewise linear reconstruction in the upstream neighboring cells are employed to interpolate (a;^ — |s) from the 
cell centered 0”''’", where the neighboring cells are identified by the direction of the particle velocity |. To demonstrate 
this procedure, here we consider a general case as illustrated in Fig. AB is a cell interface with its center locating 
a,t Xf and the unit normal vector Uf pointing from cell P to cell N. The distribution function ^'^{xf — ^s,^,tn) is 
evaluated as 


= 0^(a^c,l) + {xf -xc- |s)'i/'(a:c,l) ■ V<))+(a:c,I), (22) 

where C stands for P if | • n,/ > 0, or A otherwise. The gradient at the cell center is determined using the least 
square method. For instance, the gradient of cell P is evaluated as 


(V0+)p 




d, 





(23) 


where the tensor G is defined as 


i 


( 24 ) 
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with di being the spatial vector from P to its ith adjacent cell center Ni, and uJi = l/|(ii| being the weighting factor. 
The function tp{xc,^) in Eq. ( |22[ ) denotes the gradient limiter which is used to suppress numerical oscillations in 
regions with discontinuities, such as the shock layer in continuum regime. In this work, we adopt the Venkatakrishnan 
limiter m which is a typical one for flow computations on unstructured meshes. 

The time step in the DUGKS is determined by the Courant-Friedrichs-Lewy (CFL) condition. 


At = a 


Ax 


IC/I + I^I 


(25) 


where 0 < a < 1 is the CFL number. Aa; is the distance between the centers of two neighboring cells that share an 
interface. 

After ge tting 0 at face centers according to Eqs. (21) and([^, the original distribution functions 0 can be recovered 
from Eq. (20a). The macro variables at time tn+iji that used to evaluate the equilibrium distribution functions 0'^ 
are calculated from 0 as 


p= \m. pu= pE=^ ie-g+h)dt 


(26) 


and 


q = 


j-^9,with«=l/'c(c»9+ (>)<({. 


(27) 


Then the flux across each cell interface can be evaluated according to Eq. (14). Finally, the cell centered 0 can be 
advanced to the new time level according to Eq. (15). 

The updating procedures presented above are all based on continuous velocity space for convenience. In actual 
implementation, the continuous velocity space is discretized into a finite discrete velocity set like the DVM [3], 
and the distribution functions such as g and h are defined at these discrete velocity points as gi and hi. Proper 
quadrature rule such as the Newton-Cotes and Causs-Hermite, is used to approximate the moments. 


p = '^Wigi, pU = ^w,^tgi, pE =-'^Wt ^'ig^ + hi 

i i i 

where the Wi are the weight coefficients for the corresponding quadrature rule. 


(28) 


III. NUMERICAL EXAMPLES 

We will apply the proposed DUCKS on unstructured meshes to two internal and one external flows to demonstrate 
its performance in multiscale flow simulations. The first one is the two dimensional lid driven cavity flow at different 
flow regimes. The second one is a multiscale unsteady gas expansion problem in which the Knudsen number ranges 
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from 10“^ to 10. The last one is a supersonic rarefied gas flow with Mach number Ma=5 passing through a circular 
cylinder at Kn = 0.1 and 1. 

As the DUGKS is an explicit scheme, the simulations start from an equilibrium state based on given initial macro 
fields. For steady problems, such as the first or the last test case, the flow fields evolve into the hnal steady states, 
which is defined as the average relative change of the temperature held in two successive steps being less than 10“®, 
be.. 


£ = 




in+1 


_ 


E 'T'n 
i ^ i 


< 10 “ 


(29) 


where the summations are taken over all of the cells. 

In all of the tests the simulated gas is argon, with molecular mass m = 6.63“^®kg and molecular diameter d = 
4.17 X 10 ^°m. The viscosity of the gas is assumed to depend on temperature following a power-law, 

M = Mre/ J (30) 


where fXref is the viscosity at the reference temperature T^ef- Here we choose oj = 0.81, and the referenced viscosity 
is set to be that of a hard-sphere gas, as used in DSMC [T]. 


A. Cavity flow 

The two dimensional lid driven cavity flow is a standard benchmark problem for the validation of classical CFD 
methods in continuum regime. This problem has also been studied recently by Benziet al. HU using a parallel DSMC 
code at Knudsen numbers Kn = 10,1.0, 0.075 and was later used as an benchmark test case to validate the UGKS and 
DUGKS in a wide range of flow regimes [I3HII1I10]- To demonstrate that the DUGKS can recover the Navier-Stokes 
limit without resolving the mean free path scale in the continuum regime, we also simulate this case in full range of 
Knudsen numbers. 

The flow domain is a square cavity with length L = Im. The upper wall moves with a constant velocity Uw, while 
other walls are kept fixed. The temperature at the four walls is fixed at = 273K and is used as the referenced 
temperature. The walls are assumed to be fully diffusive and the boundary conditions are realized following the 
method presented in Refs. [laHU. The Knudsen number is defined as Kn = A/L, where A is the mean-free-path of 
the gas. Different Knudsen number can be achieved by adjusting the initial density Pref- 

Both rarefied and continuum flows are simulated. In rarefied regimes, three values of the Knudsen number, Kn = 
10,1 and 0.075, are considered. The velocity of the upper wall is set to be t/w = 50m/s, which is the same configuration 
as used in the DSMC and UGKS simulations [13 [IS]. For continuum flows, two Reynolds numbers are considered, 
i.e., Re = 400 and 1000. Here the Reynolds number is defined as Re = prefLUw/pref, and the corresponding Knudsen 
numbers are 3.7763x 10“^ and 1.5105x 10“"^, respectively. Furthermore, the Mach number Ma = Uw/-s/^RTref = 0-1, 
so that the flow is nearly incompressible and we can compare our results with the benchmark solutions mi based on 
the incompressible Navier-Stokes equations. 

In previous work m. the DUGKS with structured meshes has been employed to simulate the cavity flow at different 
flow regimes. Here we choose unstructured meshes to demonstrate the performance of the proposed method. Figure]^ 
presents the meshes used for flows with finite Knudsen numbers and continuum flows, respectively. Note that the 
mesh in Fig.j^b) is a hybrid mesh with quadrilateral cells near the walls, which performs better than a pure triangular 
mesh in capturing the boundary layer effect that is important in continuum flows. 

The discretization of velocity space and quadrature rules are chosen dependent on the Knudsen number. For 
highly rarefied flows, {i.e., Kn = 10,1), We will use the Newton-Cotes rule with 101 x 101 velocity points distributed 
uniformly in the range of [—4:^y2RT^, 4-\/2i?Tw] x [—4-\/2i?Tw, 4-\/2R7U] • For the case of Kn = 0.075, we will adopt the 
half-range Gauss-Hermit quadrature with 28 x 28 velocity points. For continuum flows, we will employ the half-range 
Gauss-Hermit quadrature rule with 16 x 16 velocity points. The CFL number is fixed at 0.8 in all simulations. 

Figures [3][^ present the temperature held, heat flux, and velocity {U, V) on the vertical and horizontal center lines, 
together with the DSMC solutions for the cases of Kn = 10,1 and 0.075, respectively. It can be seen that the present 
results agree well with DSMC results. It is interesting to note that the direction of the heat flux is inconsistent with 
the temperature gradient in each case, suggesting that the Fourier law breaks down, even at the Knudsen number as 
small as 0.075. 

Figures and show the streamline and velocity profiles for the cases of Re = 400 and 1000 respectively. The 
benchmark solutions m are also included for comparison. We can see that even though the cell size is much larger 
than the mean free path in these cases,the DUGKS results are still in very close agreement with the benchmark data. 
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FIG. 2: Meshes for the cavity flow, (a) Kn = 10,1 and 0.075. (b) Re = 400 and 1000. 




Y. X 

(c) 


FIG. 3: Results of the cavity flow at Kn = 10. (a) temperature contours, black line: DSMC, white line with background: 
DUGKS, (b) heat flux: blue solid line with arrows: DSMC, red dashed line : DUCKS, (c) U-velocity along vertical center line 
and V-velocity along horizontal central line 




Y. X 

(c) 


FIG. 4: Results of the cavity flow at Kn = 1. (a) temperature contours, black line: DSMC, white line with background: 
DUGKS, (b) heat flux: blue solid line with arrows: DSMC, red dashed line : DUGKS, (c) U-velocity along vertical center line 
and V-velocity along horizontal central line 
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FIG. 5: Results of the cavity flow at Kn = 0.075. (a) temperature contours, black line: DSMC, white line with background: 
DUGKS, (b) heat flux: blue solid line with arrows: DSMC, red dashed line : DUGKS, (c) U-velocity along vertical center line 
and V-velocity along horizontal central line 


So the DUGKS recovers the Navier-Stokes solutions in the continuum limit. We would also like to point out that for 
most traditional DVM methods, the numerical dissipation is proportional to cell size due to the splitting treatment 
of particle transport and collision processes. This may lead to significant errors for unstructured meshes as the cell 
size changes dramatically. The above results indicate that the DUGKS can avoid this difficulty with the coupled 
treatment of particle transport and collision. 



(a) 

FIG. 6: Results of the cavity flow at Re = 400, Kn = 3. 
central line and V-velocity alone horizontal central line 



(b) 


X 10 (a) Velocity streamline (b) U-velocity alone vertical 


B. Multiscale expansion flow between two connected cavities 

In the above subsection, the cavity flow at specific regimes have been simulated. Now we consider a gas expansion 
flow between two connected cavities with different initial pressures. This flow is an unsteady multiscale problem where 
different flow regimes appear in a single run. The flow configuration is sketched in Fig. Two square cavities A and 
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(a) (b) 

FIG. 7: Results of the cavity flow at Re = 1000, Kn = 1.5105 x 10“"^. (a) Velocity streamline (b) U-velocity alone vertical 
central line and V-velocity alone horizontal central line 


B connected by a channel are initially maintained at different pressure and separated by a diaphragm at the middle 
of the channel. The height of the cavity is L = Im, and the length and width of channel are L and H with H = L/8. 
Initially, the temperature of the gas in the system and that of the solid walls are maintained at 273K, which is used as 
the reference temperature. The initial Knudsen numbers at cavity A and cavity B are Kn^i = 0.001 and Kn^ = 10, 
and the corresponding pressures are Pa = 48.78Pa and Pb = 0.04878Pa, respectively. At time t = 0, the diaphragm 
is removed suddenly, and then the gas starts to expand from the left cavity to the right one. We are interested in the 
dynamic behavior of the gas during the expansion process. 

The mesh for this case is shown in Fig. As significant flow variations can take place in cavity B, the mesh 
is much finer there. While in cavity A, the flow changes slowly and the mesh is relatively coarser. Note that like 
the continuum cavity flow, the cell size in cavity A is much larger than the mean free path there. The correctness 
of using such a coarser mesh in cavity A is granted by the AP property of the DUGKS. To account for the highly 
non-equilibrium effect in cavity B at the early stage, we use a 101 x 101 mesh distributed uniformly in the range of 
[—7\/2RT^, 1\/2RT^] x [—7\/2RT^, 7\/2RT^] for the velocity space discretization, and the Newton-Cotes quadrature 
rule is used for the numerical integration. Note that a wider bound is used for the discrete velocity than that used in 
the cavity flow simulations to account for the supersonic flow behavior in the channel and cavity B in the early stage 
of expansion. In our simulations, the CFL number is set to be 0.8. 

We first define a characteristic time of the system as = L/\/2RT^, and the flow fields at different times are 
measured. The local Mach number, pressure, and streamlines at times t/tc = 1 and 4 are presented in Figures I0p3 


From Fig. 10 and Fig. 11 we can see that the shock wave just reaches the center of cavity B at time t/tc = 1- At this 
time the gas is still very rarefied there, such that the viscous effect can be neglected without vortex formation. As 
the gas flows into cavity B gradually, the pressure in cavity B grows up with time, but the pressure ratio between the 
two cavities is still high enough to form a supersonic nozzle flow in the channel, and the initial shock wave disappears 
and two symmetry vortexes appear in cavity B at a later time. 

To get a detail information of the evolution process of the expansion, we show the temperature, U-velocity and 
pressure profiles along the horizontal center line of the cavities and the channel at times t/tc = 0.013,0.1,1, 2,3 and 
4 in Fig. 14 It can be seen that at the early stage {t/tc = 0.013 and 0.1), the shock wave propagates in the channel, 
and the flow variables change sharply across the shock. With time advancing, the pressure difference between the two 
cavities decreases, and the shock becomes weaker, and the flow rate decreases gradually. 

To quantify the results, the temperature, velocity and pressure profiles along the vertical center lines of the two 
cavities at different times are presented in Fig. and Fig. respectively. Here only the results at the upper 
half { Q < y < L/2 ) of the domain are shown owing to the symmetry of the flow. From Fig. [T^b), we can see 
that a counterclockwise vortex develops in the upper half of cavity B, which enhances heat convection in the gas. 
Consequently, the temperature held becomes uniform gradually, as indicated in Fig. |15[a). 
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r„ = 273K 



Cavity A 



Cavity B 

KnA = 0.001 
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FIG. 8: Gas expansion between two cavities connected by a channel 



X 

FIG. 9: Mesh for the gas expansion case 


On the other hand, the flow in cavity A changes only slightly. The temperature and pressure are almost uniform 
at each time. With the decreasing of pressure in the cavity, the temperature reduces as the internal energy being 
converted to the kinetic energy. 



FIG. 10: Mach number contours for the gas expansion problem at time t/tc = 1 
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FIG. 11: Pressure contours and streamlines for the gas expansion problem at time tjtc. = 1 



FIG. 12: Mach number contours for the gas expansion problem at time tjta = 4 


C. Supersonic flow passing through a circular cylinder 


To further demonstrate the performance of the DUGKS on unstructured meshes for high speed non-equilibrium 
external flows, we simulate the rarefied gas flows passing through a circular cylinder. It is noted that this problem was 
also studied by Huang et al. m using the UGKS method. We here adopt the same configuration and parameters as in 
their simulations. The free-stream Mach number is Maoo = 5, and the radius of the cylinder which is r = 0.01m. Two 
Knudsen numbers are considered (Knco = Aoo/r = 0.1 and 1). The free-stream gas temperature is = 273K and 
is used as the referenced temperature. The surface of the cylinder maintains a constant temperature at = 273K, 
and full diffusive boundary condition is assumed. The outer boundary of the computational domain is a circle with a 
diameter Do = 22r, and forms a concentric annular with the surface of the cylinder. The distribution functions coming 
to the computational domain from the outer boundary are set to the equilibrium state based on the free-stream flow 
condition. 


Hybrid meshes are adopted again for this test case (see Fig. 17). Locally refined quadrilateral cells are used near 


the cylinder to resolve the boundary layer. We note that the mesh resolution in the normal direction of the cylinder 
wall should be fine enough near the cylinder to capture the large gradients correctly in the boundary layer. For the 
case of Kn = 0.1, the mesh spacing around the cylinder wall is finer than that for Kn = 1 (see Fig. [l7|(b)) since 
the boundary layer become thinner as Kn goes down. It should be pointed out that the fine resolutions around the 
cylinder wall are only used to capture the large gradients of the flow field but not to resolve the mean free path scale. 
Actually, based on the posterior estimation, the mesh spacing around the stagnation point for the case of Kn = 1 is 
about 2 times that of the mean free path there. 

In our computations, the velocity space is discretized into a set of uniform spaced 89 x 89 points in the range of 
[—15y/2RToo, 15-\/2i?Too] x [—Iby/THT^, 15y/2RToo], and the Newton-Gotes quadrature rule is used for the numerical 
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FIG. 13: Pressure contours and streamlines for the gas expansion problem at time tjtc. = 4 
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FIG. 14: Temperature (a), horizontal velocity (b), pressure (c) and Mach number (d) along the horizontal center line across 
the cavities and the channel at different times for the gas expansion problem. 
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FIG. 15: Temperature (a), horizontal velocity (b), vertical velocity (c) and pressure (d) along the vertical center line (upper 
half) of the cavity B at different times for the gas expansion problem. 


integration. To validate our simulation results, we use the open source dsmcFoam solver [52] to obtain the DSMC 
results under the same flow conditions and computational domain. For the case of Kn = 1, the total number of DSMC 
particles is about 0.58 million and the time step is 1 x 10“^s. For the case of Kn = 0.1, the total number of DSMC 
particles is about 1.11 million and the time step is 4 x 10“®s. 

The contours of temperature and Mach number for the case of Kn = 0.1 are shown in Fig. also included are 
the DSMC solutions. The temperature and U-velocity prohle along the stagnation line are shown in Fig. Clearly 
we can see that both the temperature and Mach number distributions of the DUCKS results agree with those of the 
DSMC results perfectly. Flowever, there are some discrepancies in the front of the bow shock, which can be seen 
more clearly in the temperature profile. This is due to the intrinsic defect of the Shakhov model used in the current 
DUCKS [Ildl], where the collision frequency is independent of particle velocities. Despite of the small deviations, 
the temperature agrees well with the DSMC results downstream the shock, and the heat flux, normal pressure and 
shear stress distribution along the cylinder’s surface predicted by the DUCKS agree with the DSMC results quite 
well, as shown in Fig. 

For the case of Kn = 1, the temperature and Mach number distributions are presented in Fig. 21 The temperature, 
U-velocity and density profile along the stagnation line are shown in Fig. 


22 These results show that the DUCKS 


results agree with the DSMC results as well. The discrepancies in the front of the bow shock are slightly more obvious. 
This is because with the increasing of Kn, the non-equilibrium effects get stronger, thus the Shakhov model deviates 
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FIG. 16: Temperature (a) and pressure (b) along the vertical center line (upper half) of the cavity A at different times for the 
gas expansion problem. 



(a) (b) 

FIG. 17: Meshes for the flow past a cylinder, (a) Global view of the mesh, (b) local view of the meshes around the cylinder 
surface, upper: Kn = 1, lower: Kn = 0.1. 


more from the full Boltzmann collision kernel. However, the heat flux, normal pressure and shear stress along surface 
of cylinder predicted by DUGKS are still quite satisfactory in comparison with the DSMC results as shown in Fig. 
These results demonstrate that although the Shakhov model has some intrinsic defects, the DUGKS based on it can 
still give rather satisfactory predictions, particularly the flow behaviors near the body. The DUGKS can be a very 
useful engineering tool for hypersonic rarefied flow applications, especially in the regime Kn <0.1. 
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(a) (b) 

FIG. 18: Temperature (a) and Mach number (b) distribution for the flow past a cylinder at Kn=0.1. Solid white line with 
colored background: DUGKS, dashed black line: DSMG 


IV. CONCLUDING REMARKS 

In this paper, the DUGKS based on the Shakhov model developed recently [16] is extended to unstructured meshes. 
The key feature of DUGKS is that the discrete characteristic solution of the kinetic equation is used in the modeling 
of the distribution function at a cell interface. Due to the coupled treatment of the particle collision and transport 
processes, the method has the asymptotic preserving (AP) properties for the capturing NS solutions in the continuum 
flow regime. Linear reconstruction and gradient limiter are employed in the reconstruction to attain the second-order 
accuracy. 

The performance of the DUGKS on unstructured meshes has been explored by several examples covering different 
flow regimes from low speed microflows to hypersonic rarefied flows. In the transitional and slip regimes, good 
agreements between the results of current scheme and the DSMC solutions are observed; In the continuum regime, 
the DUGKS results obtained on a coarse mesh without resolving the mean free path agree with the benchmark solution 
based on the Navier-Stokes equations very well. Thus the AP property of the DUGKS for the Navier-Stokes limit is 
demonstrated. The merit of this property is important for a multicale unsteady gas expansion problem that involves 
both continuum and rarefied regions. As the mesh size in the continuum region can be much larger than the mean free 
path scale, the overall computational cost for DUGKS can be largely reduced in comparison with the DSMG method 
and the traditional DVM. Since the DUGKS is a direct modeling multiscale method [53] method, with the mesh size 
and time step being a few times of particle mean free path and collision, the physical solutions are not sensitive to 
individual particle collision. The DUGKS based on kinetic model equation can be faithfully used in real engineering 
applications. 
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FIG. 19: Temperature (a), velocity (b) and density (c) profiles alone the stagnation line for flow past a cylinder at Kn=0.1. 
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FIG. 20: Heat flux (a), pressure (b) and shear stress (c) alone the surface for the flow past a cylinder at Kn=0.1. 
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(a) (b) 

FIG. 21: Temperature and Mach number distributions for the flow past a cylinder at Kn=l. Solid white line with colored 
background: DUGKS, dashed black line: DSMG 
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FIG. 22: Temperature (a), velocity (b) and density (c) profiles alone the stagnation line for flow past a cylinder at Kn=l. 
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FIG. 23: Heat flux (a), pressure (b) and shear stress (c) alone the surface for the flow past a cylinder at Kn=l. 








21 


References 


[1] G.A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Glarendon, Oxford, 1994. 

[2] J. Fan, G. Shen, Statistical Simulation of Low-Speed Rarefied Gas Flows, J. Comput. Phys., 167 (2001) 393-412. 

[3] T.M. Homolle, N.G. Hadjiconstantinou, A low-variance deviational simulation Monte Carlo for the Boltzmann equation, 
J. Comput. Phys., 226 (2007) 2341-2358. 

[4] J.Y. Yang, J.C. Huang, Rarefied Flow Computations Using Nonlinear Model Boltzmann Equations, J. Comput. Phys., 120 
(1995) 323-339. 

[5] L. Mieussens, Discrete-velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric 
geometries, J. Comput. Phys., 162 (2000) 429-466. 

[6] Z.H. Li, H.X. Zhang, Gas-kinetic numerical studies of three-dimensional complex flows on spacecraft re-entry, J. Comput. 
Phys., 228 (2009) 1116-1138. 

[7] E.F. Toro, Riemann solvers and numerical methods for fluid dynamics. Springer, Heidelberg, 2009. 

[8] M. Bennoune, M. Lemon, L. Mieussens, Uniformly stable numerical schemes for the Boltzmann equation preserving the 
compressible NavierStokes asymptotics, J. Comput. Phys., 227 (2008) 3781-3803. 

[9] F. Filbet, S. Jin, A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources, J. 
Comput. Phys., 229 (2010) 7625-7648. 

[10] K. Xu, J.C. Huang, A unified gas-kinetic scheme for continuum and rarefied flows, J. Comput. Phys., 229 (2010) 7747-7764. 

[11] K. Xu, J.C. Huang, An improved unified gas-kinetic scheme and the study of shock structures, IMA J. Appl. Math., 76 
(2011) 698-711. 

[12] J.C. Huang, K. Xu, P.B. Yu, A Unified Gas-Kinetic Scheme for Continuum and Rarefied Flows H: Multi-Dimensional 
Cases, Commun. Comput. Phys., 12 (2012) 662-690. 

[13] J.C. Huang, K. Xu, P.B. Yu, A Unified Gas-Kinetic Scheme for Gontinuum and Rarefied Flows HI: Microflow Simulations, 
Gommun. Gomput. Phys., 14 (2013) 1147-1173. 

[14] S.Z. Ghen, K. Xu, A comparative study of an asymptotic preserving scheme and unified gas-kinetic scheme in continuum 
flow limit, J. Comput. Phys., 288 (2015) 52-65. 

[15] Z.L. Guo, K. Xu, R.J. Wang, Discrete unified gas kinetic scheme for all Knudsen number flows: Low-speed isothermal 
case, Phys. Rev. E, 88 (2013) 033305. 

[16] Z.L. Guo, R.J. Wang, K. Xu, Discrete unified gas kinetic scheme for all Knudsen number flows: 11. Gompressible case, 
arXiv:1406.5668 (2014). 

[17] E. Shakhov, Generalization of the Krook kinetic relaxation equation, Fluid Dynamics, 3 (1968) 95-96. 

[18] V. Venkatakrishnan, Gonvergence to Steady State Solutions of the Euler Equations on Unstructured Grids with Limiters, 
J. Comput. Phys., 118 (1995) 120-130. 

[19] B. John, X.J. Gu, D.R. Emerson, Effects of incomplete surface accommodation on non-equilibrium heat transfer in cavity 
flow: A parallel DSMG study, Comput. Fluids, 45 (2011) 197-201. 

[20] P. Wang, L.H. Zhu, Z.L. Guo, K. Xu, A comparative study of LBE and DUGKS methods for nearly incompressible flows, 
Gommun. Gomput. Phys., 17 (2015) 657-681. 

[21] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid 
method, J. Comput. Phys., 48 (1982) 387-411. 

[22] T.J. Scanlon, E. Roohi, C. White, M. Darbandi, J.M. Reese, An open source, parallel DSMC code for rarefied gas flows in 
arbitrary geometries, Comput. Fluids, 39 (2010) 2078-2089. 

[23] K. Xu, Direct Modeling for Computational Fluid Dynamics: Construction and Application of Unified Gas-Kinetic Schemes, 
World Scientific Publishing, 2015. 



